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INTRODUCTION 


Propellant  formulators  and  ballisticians  have  been  using  thermodynamic  codes 
such  as  Blake,  Nasa-Lewis,  and  more  recently  MCVEC  (refs  1  through  3,  respectively) 
to  predict  various  thermodynamic  parameters  and  combustion  products  for  gun  firings. 
The  critical  assumption  used  in  applying  these  codes  is  that  sufficient  time  has  elapsed 
during  the  interior  ballistic  cycle  to  achieve  a  state  of  minimum  free  energy  (chemical 
equilibrium)  within  the  combustion  chamber.  The  validity  of  this  assumption,  however, 
has  not  been  determined,  and  as  a  result,  should  be  thoroughly  investigated.  In  theory, 
this  could  be  done  by  analyzing  the  chamber  combustion  products  and  comparing  them 
with  the  equilibrium  code  data. 

Since  experimental  gun-shot  data  are  not  always  available,  theoretical  methods 
involving  chemical  kinetics  can  be  employed  to  map  the  reaction  paths  and  predict  the 
species  generated  during  the  interior  ballistic  cycle.  The  magnitude  of  this  undertaking, 
however,  requires  the  use  of  computers  to  perform  the  calculations  efficiently.  Several 
codes  which  model  chemical  kinetics  are  available  for  this  purpose,  but  unfortunately 
these  existing  computer  codes  are  written  for  either  main-frame  computer  application  or 
require  long  run  times. 

Recent  developments  in  using  virtual  memory  on  personal  computers  have  guided 
the  initial  intention  of  this  effort  to  develop  a  kinetic  code  which  will  run  on  a  desk-top 
computer  and  yet  not  require  extensive  running  time.  Using  such  a  system,  it  may  be 
possible  to  perform  an  approximate  gun  calculation  on  a  personal  computer.  The 
ultimate  goal,  however,  is  to  model  the  propellant  combustion  including  both  the  reac¬ 
tion  kinetics  and  propellant  hydrodynamics.  Therefore,  a  chemical  reaction  kinetics 
code  was  written  for  use  on  a  personal  computer  which  would  eventually  form  part  of 
the  larger  program.  In  addition,  such  a  code  would  be  useful  as  a  tool  to  supplement 
the  chemical  equilibrium  codes  currently  in  use. 

To  test  the  validity  of  the  assumptions  incorporated  into  the  algorithm  of  the  code, 
the  oxidation  of  hydrogen  at  various  conditions  was  run  and  compared  against  the 
results  from  the  equilibrium  code.  It  is  difficult  to  undertake  a  comparison  with  steady 
state  experimental  results  as  most  experiments  have  been  performed  at  between  1000 
K  and  2000  K  where  the  main  product  is  water. 

The  initial  propellant  selected  to  develop  this  kinetic  code  was  the  HAN  based  liquid 
propellant  system  (LP1845  and  LP1846).  This  selection  was  based  on  the  current 
interest  in  HAN  based  regenerative  liquid  propellant  guns  and  the  apparently  simple 
molecules  involved  in  LP  formulations:  HAN  (NH„OHNO„),  TEAN  (C,H,0),NHNO„,  and 

water  (H^O). 
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Only  the  oxidation  of  anhydrous  HAN  is  discussed  in  this  report.  The  kinetic  code 
results  are  compared  to  equilibrium  code  data  as  a  first  approximation  goodnesss-of-fit. 
However,  since  the  test  of  any  modeling  scheme  should  be  agreement  with  experimen¬ 
tal  data,  the  kinetic  code  results  are  also  compared  with  experimental  thermal  decom 
position  data. 


METHOD  OF  CALCULATION 

The  equations  required  for  the  chemical  kinetics  of  the  fluid  system  will  involve 
variable  temperature  and  density.  The  present  calculations  were  undertaken  at  con¬ 
stant  density  as  this  is  most  easily  compared  with  the  minimum  free  energy  calculations. 
The  changes  required  to  allow  variable  density  are  minimal.  The  basic  equation  for  the 
change  of  the  concentration  c  of  the  ith  species  is 

§=ip,r,  (1) 

where  (3.^  is  the  ith  species  in  the  jth  equation  and  r.  is  defined  as  the  net  rate  by  the 
mass-action  equation 

rj  =  kf  nci  -  kr  nci  (2) 

where  the  forward  reaction  of  rate  k,  is  related  to  the  reverse  reaction  rate  k  by  the 
equilibrium  constant 

^  =  Ke,  (3) 

The  equilibrium  constant  used  expresses  the  ratio  between  the  number  of  moles  of 
products  and  the  number  of  moles  of  reactants.  In  the  case  where  the  number  of  moles 
of  products  and  reactants  are  the  same,  then  the  usual  equilibrium  constant.  In 

the  case  that 


AN  =  IN  -  IN  (4) 

p  r 

is  not  zero,  then 

K.,  =  ^(RT) (5) 
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In  the  case  studied  here,  the  pressure  P  will  rise  according  to  the  relation  for  ideal 
gases 


dP 

dt 


RTI 


dCi 

dt 


(6) 


Equation  1  is  used  to  integrate  the  change  in  concentration  in  time.  Conceptually,  the 
integration  ot  this  derivative  Is  quite  straightforward,  but  the  real  difficulties  are  imposed 
by  the  limitation  in  accuracy  of  the  computer.  The  first  problem  is  the  fact  that  the 
integration  .scheme  may  become  unstable  with  time  steps  which  are  too  large.  This  can 
be  overcome  by  various  methods.  The  most  commo  ,  method  involves  calculating  the 
high  order  derivatives  (for  example,  the  first  four  derivatives  and  enforcing  continuity  in 
all  derivatives).  A  major  problem  with  this  method  is  that  sudden  changes  cannot  be 
made  in  a  variable.  These  changes  must  be  undertaken  gradually.  This  is  satisfactory 
on  a  main-frame  computer  where  calculation  time  is  not  the  chief  concern.  When 
solving  the  problem  on  a  PC  with  reasonable  run  time,  it  is  necessary  to  be  able  to 
make  approximations  which  require  sudden  changes  in  the  function  in  order  to  con¬ 
serve  certain  quantities  (e.g.,  atom  numbers).  Therefore,  a  semi-implicit  time  advanced 
scheme  (ref  4)  was  chosen. 

If  equation  is  written  as 


^  -  P 

dt 


Xi  L| 


(7) 


where  P^  are  the  production  terms  and  L.  are  the  loss  terms  for  the  species  x,.  The 
usual  finite  difference  form  of  the  semi-implicit  scheme  is 


Xj  (j  +  1 )  -  Xj  (j) 
t 


=  Pi  (j)  -  Xi  (j  1  )Li  (j) 


(8) 


where  the  right  hand  side  can  be  written  as 


1 

1  +  tLi  (j) 


f(j).  Now  equation  7  can  be 


approximated  by 


dx 


Sf  =Pi(t,)-Xi(t)U(tj) 


(9) 


integrating  equation  9  we  obtain 

X|  (j  +  1 )  -  Xj  (j) 

t 


1  -  e 


IL: 


tLi  (j) 


t(j) 


(10) 
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where  t  is  the  time  step.  This  scheme  is  very  stable  and  the  original  paper  should  be 
consulted  for  details. 

The  second  problem  involves  dealing  with  a  large  range  in  the  time  derivatives.  In 
this  case,  the  range  is  quite  large  because  the  reaction  rates  vary  from  10E-18  to 
10E+20  sec  V  These  equations  are  known  as  stiff,  and  various  schemes  have  been 
proposed  to  handle  them.  For  example,  if  the  computer  can  retain  12  significant  figures 
and  the  time  step  is  chosen  as  10E-8  sec,  then  only  the  last  couple  of  digits  of  the 
smallest  derivative  will  be  involved  while  the  largest  derivative  will  not  be  involved  The 
number  of  time  steps  to  perform  an  integration  will  be  prodigious  for  a  PC. 

The  critical  mechanism  for  propellant  combustion  is  to  decompose  the  initial  in¬ 
gredients  into  reactive  gaseous  species  that  will  drive  the  reactions.  This  is  the  overall 
rate-determining  step  and  the  sequence  of  subsequent  time  steps  must  be  chosen 
around  this  rate.  This  is  necessary  in  order  to  make  the  propellant  decompose  to  start 
the  reaction. 

As  a  first  approximation,  the  initial  shortest  time  step  for  the  reaction  (10E-18  sec)  is 
arbitrarily  set  to  equal  one  tenth  the  value  of  the  fastest  reaction  characteristic  time  (in 
this  case  10E-17  sec).  A  characteristic  time  for  a  given  reaction  is  defined  as  the 
maximum  time  required  to  use  up  one  of  the  reactants  and  is  approximately  the  recipro¬ 
cal  of  reaction  rate.  The  next  time  step  is  determined  as  a  function  of  the  propellant 
decomposition  equations.  Any  equation  involving  the  propellant  decomposition  may  be 
used  for  this  although  some  experimentation  is  required  to  find  the  most  efficient  equa¬ 
tion.  The  fastest  time  step  is  then  repeated  and  finally  the  last  time  step  is  a  function  of 
the  propellant  decomposition  equations.  The  program  is  flexible  so  that  the  time  steps 
can  be  altered  during  a  run  if  desired  to  allow  for  special  cases. 

The  technique  used  here  is  related  to  that  adopted  in  reference  5  although  the 
implementation  is  different.  A  variety  of  time  steps  are  chosen  so  that  all  derivatives 
may  contribute  in  varying  degrees.  For  example  10  time  steps  of  length  10E-17  sec 
then  30  of  length  10E-13  sec,  followed  by  10  at  10E-17  sec  and  3  at  10E-12  sec,  10  at 
10E-17  sec  and  2  at  10E-1 1  sec.  The  sequence  is  then  repeated.  The  time  steps  are 
assigned  to  each  reaction  as  a  function  of  the  characteristic  time  for  any  component  in 
the  reaction.  Therefore,  in  starting  a  run,  it  is  necessary  to  run  a  few  time  steps  with 
very  short  times  to  obtain  an  approximation  of  the  amounts  of  various  gaseous  sub¬ 
stances  formed.  The  time  steps  are  revised  at  frequent  intervals  usually  every  three 
times  the  number  of  steps  in  each  cycle.  For  example,  in  the  above  case,  it  would  be 
every  195  time  steps. 
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Finally  it  was  found  that  atom  conservation  becomes  difficult  when  using  a  PC. 
This  was  overcome  by  performing  a  count  of  atoms  at  the  end  of  each  time  step  and 
adjusting  the  number  of  free  carbon,  hydrogen,  oxygen,  and  nitroger.  atoms  to  conserve 
atoms.  These  were  chosen  since  they  are  very  reactive  and  quickly  adjust  in  the  follow¬ 
ing  time  step.  The  proposed  scheme  appears  to  work  well  for  propellants  due  to  the 
relative  simplicity  of  the  reactions  involved  and  the  limited  number  of  different  species 
involved. 


RESULTS 

The  initial  tests  were  undertaken  on  a  mixture  of  one  mole  of  and  half  mole  of 
O2.  The  program  was  tested  against  results  obtained  from  the  chemical  equilibrium 

code  MCVEC.  The  results  were  obtained  at  1000  K/6  atm  and  at  the  adiabatic  flame 
temperature  of  3600  K/20  atm.  The  results  are  shown  in  table  1 . 


Table  1 .  Comparison  of  minimum  free  energy  and  reaction  kinetics  for 
H2  -I-  2O2  at  1000  K  and  3600  K 

Mole  fractions 


1000  K  3600K 


Species 

MCVEC  code 

Kinetic  code 

MCVEC  code 

Kinetic  code 

H,0 

1 .0000 

0.9907 

0.68 

0.73 

<10E-5 

0.0015 

0.08 

0.16 

H 

<10E-5 

0.002 

0.05 

0.02 

OH 

<10E-5 

0.004 

0.20 

0.002 

0 

<10E-5 

0.0004 

<10E-5 

0.0001 

0, 

<10E-5 

0.0014 

<10E-5 

0.007 

At  1000  K,  the  kinetic  and  equilibrium  codes  are  in  good  agreement  (table  1),  but  a 
noticeable  difference  between  the  two  codes  is  observed  at  the  higher  temperature 
where  the  equilibrium  code  yields  less  H^O  and  more  OH  than  the  kinetic  code.  Since 
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published  reaction  rates  for  the  water-OH  reactions  vary  over  a  wide  range,  a  mean 
value  was  chosen  for  the  present  test.  The  time  variation  of  the  species  at  3600  K  is 
shown  in  figure  1 . 

The  program  was  then  tested  using  the  experimental  results  (ref  6)  for  HAN.  The 
HAN  was  heated  from  300  K  to  650  K  in  10  seconds.  The  kinetic  code  allows  for  this 
change  of  temperature  as  a  linear  variation.  The  implementation  of  the  code  requires 
the  use  of  various  simple  reactions  as  well  as  an  initial  breakdown  scheme  for  the  HAN. 
The  initial  breakdown  equation  used  was 

=  jN^O  +  H^O  +  O  +  ^HNOg  +  JnH^  (11) 

The  final  equilibrium  results  were  obtained  using  the  MCVEC  code  for  ideal  gases. 
The  proposed  reaction  kinetic  scheme  was  used  to  model  the  heating,  and  the  results 
are  shown  in  figure  2.  A  comparison  of  the  final  states  using  the  three  methods  is 
shown  in  table  2.  A  comparison  of  the  final  states  using  the  three  methods  is  shown  in 
table  2.  The  experimental  results  are  in  volume  fractions,  and  the  calculated  ones  are 
mole  fractions.  In  addition,  the  experimental  values  do  not  include  species  which  are  IR 
inactive  nor  was  the  water  volume  measured.  Therefore,  the  experimental  values  only 
show  the  relative  proportions  of  the  different  species. 


Table  2.  Comparison  of  experimental,  equilibrium,  and  kinetic  results  for  HAN 


heated  to  650  K  and  1  atm 

Experimental 

Mole  fractions 

Species 

volume  fraction 

MCVEC  code 

Reaction  kinetic 

HNO3 

0.76 

<10E-5 

0.135 

NO2 

0.2 

<10E-5 

<10E-02 

N,0 

0.02 

<10E-05 

<10E-02 

N, 

0.250 

0.2236 

H3O 

0.500 

0.5106 

O2 

0.250 

0.125 

NO 

<10E-05 

0.0011 
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The  experimental  results  shown  in  table  2  considerable  HNO3  compared  to  the  NO^ 

formed.  This  is  probably  due  to  the  way  the  HAN  decomposed  which  strongly  suggests 
that  the  equilibrium  calculations  are  very  long-time  results;  much  longer  times  than  are 
of  interest  in  gun  calculations.  The  amount  of  HNO3  shown  by  the  kinetic  code  is  not 

necessarily  incorrect  since  the  amount  of  HNO3  from  the  experimental  results  can  only 

be  compared  to  the  nitrogen  oxides  formed.  As  these  were  not  produced  by  the  kinetic 
code,  a  comparison  cannot  be  made.  It  suggests  that  considerable  work  is  needed  to 
determine  correct  reaction  rates  for  HNO3  because  only  limited  data  are  available.  The 

amount  of  water  calculated  by  both  codes  is  approximately  the  same.  The  kinetic  code 
predicts  more  and  less  NO2  than  the  equilibrium  code.  It  appears  that  measure¬ 
ments  of  the  water  formed  would  be  a  useful  measure  to  judge  the  accuracy  of  the 
calculations.  The  decomposition  products  of  HAN  when  heated  from  300  K  to  1600  K 
and  15  atms  were  calculated  and  measured  by  all  three  methods.  The  results  are 
shown  in  table  3. 


Table  3.  Comparison  of  experimental,  equilibrium,  and  reaction  kinetic  results 
for  HAN  heated  to  1600  K  and  15  atm 


Experimental 

Mole  fractions 

Species 

volume  fraction 

MCVEC  code 

Reaction  kinetic 

HNO3 

0.2 

<10E-05 

0.114 

NO2 

0.28 

<10E-05 

<0.01 

N^O 

0.5 

<10E-05 

<0.01 

N. 

0.249 

0.22 

H2O 

. 

0.499 

0.525 

0, 

0.249 

0.05 

NO 

0.0015 

0.0107 

Again  the  equilibrium  calculations  failed  to  show  the  formation  of  HNO3  while  the 

kinetic  code  produced  a  significant  amount.  The  water  produced  in  both  the  calcula¬ 
tions  was  approximately  the  same.  Problems  with  the  formation  of  nitrogen  oxides  were 
apparent  in  the  kinetic  code. 
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CONCLUSIONS 


The  kinetic  code  appears  to  clearly  show  that  it  agrees  well  with  the  equilibrium 
code  for  long  run  times.  There  remains  some  doubt  about  the  formation  of  HNO3  and 

nitrogen  oxides  when  compared  to  experimental  results.  This  suggests  that  the  equilib¬ 
rium  codes  are  not  always  suitable  for  gun  calculations  because  it  may  require  times 
much  longer  than  milliseconds  for  the  gun  reactions  to  reach  equilibrium. 


RECOMMENDATIONS 

Further  development  of  this  code  will  involve  studies  of  the  wet  HAN,  dry  TEAM,  wet 
TEAN,  dry  HAN/TEAN,  and  LP1 845/1 846.  There  still  remains  additional  research  to  be 
undertaken  on  optimizing  the  code  and  allowing  for  certain  reaction  paths  to  be  blocked 
due  to  three  body  effects.  The  inclusion  of  transport  properties  will  be  necessary  before 
the  hydrodynamic  model  can  be  developed. 
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